/*
 * Measure.cpp
 *
 *  Created on: Mar 12, 2012
 *      Author: snirgaz
 */

#include "Measure.h"

#define PI 3.14159265

// Winding Measure
Measure::Measure() {

}
Measure::~Measure() {

}

WindingEnergyMeasure::WindingEnergyMeasure(SSEConf *sseConf,
		QmcParameters &qmcParameters) {
	sseConf_ = sseConf;
	ScalarObs<double>::InitData initDataScalar = {
			qmcParameters.numOfBins.getVal(), qmcParameters.binSize.getVal() };
	windX_.init(initDataScalar);
	windY_.init(initDataScalar);
	energy_.init(initDataScalar);
	cv_.init(initDataScalar);
}
void WindingEnergyMeasure::DoMeasurement() {
	boost::array<int, 2> wind;
	wind.fill(0);
	boost::for_each(sseConf_->opList, [&](Op & op) {
		//for(auto it_o=sseConf_->opList.begin();it_o!=sseConf_->opList.end();it_o++)
		//	Op& op=*it_o;
			BondToSite bondToSite;
			if (op.type == OFF) {
				bondToSite = sseConf_->bondsArray[op.bond];
				int d = (op.bond < sseConf_->latSize) ? 0 : 1;
				if (sseConf_->spin[bondToSite[0]] == -1)
				wind[d]++;
				else
				wind[d]--;
				boost::for_each(bondToSite, [&](int &site) {
							sseConf_->spin[site]=-sseConf_->spin[site];
						}
				);
			}
		});
	double windsqrX = 0.5 * (wind[0] * wind[0]);
	double windsqrY = 0.5 * (wind[1] * wind[1]);
	windX_ << windsqrX;
	windY_ << windsqrY;
	energy_ << double(sseConf_->expRank);
	cv_ << double(sseConf_->expRank * sseConf_->expRank);
}

bool WindingEnergyMeasure::checkDone() const {
	return windX_.checkDone();
}
void WindingEnergyMeasure::dumpToFile(hdf5RW &hdf5File) const {
	windX_.dumpToFile(hdf5File, "WINDX");
	windY_.dumpToFile(hdf5File, "WINDY");
	energy_.dumpToFile(hdf5File, "ENERGY");
	cv_.dumpToFile(hdf5File, "CV");
}
void WindingEnergyMeasure::readFromFile(hdf5RW &hdf5File) {
	windX_.readFromFile(hdf5File, "WINDX");
	windY_.readFromFile(hdf5File, "WINDY");
	energy_.readFromFile(hdf5File, "ENERGY");
	cv_.readFromFile(hdf5File, "CV");
}
int WindingEnergyMeasure::getNumOfMeasure() const {
	return windX_.getNumOfMeasure();
}

// Steggered Magnetization

StaggeredMeasure::StaggeredMeasure(SSEConf *sseConf,
		QmcParameters &qmcParameters) {
	sseConf_ = sseConf;
	ScalarObs<double>::InitData initDataScalar = {
			qmcParameters.numOfBins.getVal(), qmcParameters.binSize.getVal() };
	mz_.init(initDataScalar);
	L_ = sseConf_->L;
}
void StaggeredMeasure::DoMeasurement() {
	double m = 0, msqr = 0;
	int ix = 0, iy = 0;
	for_each(sseConf_->spin.begin(), sseConf_->spin.end(),
			[&](SSEConf::SpinType spin) {
				int sign=((ix+iy)%2) ? 1 : -1;
				m+=sign*spin;
				ix=(ix+1);
				if (ix == L_) {
					ix =0;
					iy=iy+1;
				}
			});
	boost::for_each(sseConf_->opList, [&](Op & op) {
		//for(auto it_o=sseConf_->opList.begin();it_o!=sseConf_->opList.end();it_o++)
		//	Op& op=*it_o;
			BondToSite bondToSite;
			if (op.type == OFF) {
				bondToSite = sseConf_->bondsArray[op.bond];
				boost::for_each(bondToSite, [&](int &site) {
							sseConf_->spin[site]=-sseConf_->spin[site];
						});
				iy=bondToSite[0] / L_;
				ix= bondToSite[0] % L_;
				int sign=((ix+iy)%2) ? 1 : -1;
				m=m+2*sign*sseConf_->spin[bondToSite[0]];
			}
			if (op.type != EMPTY) msqr=msqr+m*m;
		});
	mz_ << 0.25 * double(msqr) / double(sseConf_->latSize * sseConf_->expRank);
}

bool StaggeredMeasure::checkDone() const {
	return mz_.checkDone();
}
void StaggeredMeasure::dumpToFile(hdf5RW &hdf5File) const {
	mz_.dumpToFile(hdf5File, "STAGGERED");
}
void StaggeredMeasure::readFromFile(hdf5RW &hdf5File) {
	mz_.readFromFile(hdf5File, "STAGGERED");
}
int StaggeredMeasure::getNumOfMeasure() const {
	return mz_.getNumOfMeasure();
}
// Lodoun Florie
HcorrMeasure::HcorrMeasure(SSEConf *sseConf, QmcParameters &qmcParameters,
		MTGen *gen) {
	sseConf_ = sseConf;
	VectorObs<double>::InitData initDataVector =
			{ qmcParameters.numOfBins.getVal(),
					qmcParameters.numOfOmegaMat.getVal(),
					qmcParameters.binSize.getVal() };
	X_.init(gen, UniRealParams(0., 1.));
	gOmega_ = ublas::zero_vector<double>(qmcParameters.numOfOmegaMat.getVal());
	HCorr_.init(initDataVector);
	numOfOmega_ = qmcParameters.numOfOmegaMat.getVal();
	numOfSites_ = sseConf_->latSize;
	beta_ = qmcParameters.beta.getVal();
}
void HcorrMeasure::DoMeasurement() {
	int expRank = sseConf_->expRank;
	Xp1_.resize(expRank);
	Xp2_.resize(expRank);
	XCorr_.resize(expRank - 1);
	tau_.resize(expRank - 1);
	std::fill(XCorr_.begin(), XCorr_.end(), 0.);
	fill(gOmega_.begin(), gOmega_.end(), 0.);
	vector<double>::iterator XPos1 = Xp1_.begin();
	vector<double>::iterator XPos2 = Xp2_.begin();
	BondToSite bondToSite;
	// Generate X
	boost::for_each(sseConf_->opList, [&](Op & op) {
		//for(auto it_o=sseConf_->opList.begin();it_o!=sseConf_->opList.end();it_o++)
		//	Op& op=*it_o;
		// Anti-parallel X -Anti-parallel Y
			if (op.type == OFF) {
				bondToSite = sseConf_->bondsArray[op.bond];
				boost::for_each(bondToSite, [&](int &site) {
							sseConf_->spin[site]=-sseConf_->spin[site];
						});
				if ((op.bond<numOfSites_) && (op.bond%2))
				*XPos2++=1./sseConf_->g;
				else
				*XPos2++=0.;
				*XPos1++=0.;
			}
			if (op.type == DIAG) {
				int strongbond = 0, pos = 0,count=0;
				for (auto itb = sseConf_->bondsArray.begin();
						itb != sseConf_->bondsArray.end(); itb++) {
					if (sseConf_->spin[(*itb)[0]] != sseConf_->spin[(*itb)[1]]) {
						strongbond += ((pos < numOfSites_) && (pos %2)) ? 1 : 0;
						count++;
					}
					pos++;
				}
				if ((op.bond%2) && (op.bond<numOfSites_))
				*XPos1++=(double(strongbond)/(double(count)))/sseConf_->g;
				else
				*XPos1++=double(strongbond)/(double(count));
				*XPos2++=0.;
			}
		});
	// Corr
	for (int p = 0; p < (expRank - 1); p++) {
		for (int m = 0; m < expRank; m++) {
			XCorr_[p] += (3. / 2.)
					* (2. * Xp1_[m] * Xp1_[(m + p + 1) % expRank]
							+ Xp1_[m] * Xp2_[(m + p + 1) % expRank]
							+ Xp2_[m] * Xp1_[(m + p + 1) % expRank]);
		}
	}
	for (vector<double>::iterator it = tau_.begin(); it != tau_.end(); it++)
		*it = X_.draw();
	std::sort(tau_.begin(), tau_.end());
	vector<double>::iterator it_tau = tau_.begin();
	vector<double>::iterator it_XCorr = XCorr_.begin();
	for (int k = 0; k < numOfOmega_; k++) {
		double temp = 0;
		vector<double>::iterator it_tau = tau_.begin();
		vector<double>::iterator it_XCorr = XCorr_.begin();
		for (; it_XCorr != XCorr_.end(); it_XCorr++, it_tau++)
			temp += (*it_XCorr) * cos(2. * PI * (*it_tau) * double(k));
		gOmega_[k] = temp;
	}
	HCorr_ << (1. / (pow(beta_, 2))) * gOmega_;
}

bool HcorrMeasure::checkDone() const {
	return HCorr_.checkDone();
}
void HcorrMeasure::dumpToFile(hdf5RW &hdf5File) const {
	HCorr_.dumpToFile(hdf5File, "HCORR");
}
void HcorrMeasure::readFromFile(hdf5RW &hdf5File) {
	HCorr_.readFromFile(hdf5File, "HCORR");
}
int HcorrMeasure::getNumOfMeasure() const {
	return HCorr_.getNumOfMeasure();
}

//

LodounFlourieMeasure::LodounFlourieMeasure(SSEConf *sseConf,
		QmcParameters &qmcParameters, MTGen *gen) {
	sseConf_ = sseConf;
	VectorObs<double>::InitData initDataVector =
			{ qmcParameters.numOfBins.getVal(),
					qmcParameters.numOfOmegaMat.getVal(),
					qmcParameters.binSize.getVal() };
	X_.init(gen, UniRealParams(0., 1.));
	gOmega_ = ublas::zero_vector<double>(qmcParameters.numOfOmegaMat.getVal());
	LF_.init(initDataVector);
	numOfOmega_ = qmcParameters.numOfOmegaMat.getVal();
	numOfSites_ = sseConf_->latSize;
	beta_ = qmcParameters.beta.getVal();
}

void LodounFlourieMeasure::DoMeasurement() {
	int expRank = sseConf_->expRank;
	Xp1_.resize(expRank);
	Xp2_.resize(expRank);
	XCorr_.resize(expRank - 1);
	tau_.resize(expRank - 1);
	std::fill(XCorr_.begin(), XCorr_.end(), 0.);
	fill(gOmega_.begin(), gOmega_.end(), 0.);
	vector<double>::iterator XPos1 = Xp1_.begin();
	vector<double>::iterator XPos2 = Xp2_.begin();
	BondToSite bondToSite;

	// Generate X
	boost::for_each(sseConf_->opList, [&](Op & op) {
		//for(auto it_o=sseConf_->opList.begin();it_o!=sseConf_->opList.end();it_o++)
		//	Op& op=*it_o;
		// Anti-parallel X -Anti-parallel Y
			if (op.type == OFF) {
				bondToSite = sseConf_->bondsArray[op.bond];
				boost::for_each(bondToSite, [&](int &site) {
							sseConf_->spin[site]=-sseConf_->spin[site];
						});
				if (op.bond>=numOfSites_) {
					*XPos2++=-1.;
				}
				else {
					if (op.bond%2) {
						*XPos2++=1./(sseConf_->g);
					}
					else
					*XPos2++=1.;
				}
				*XPos1++=0.;
			}
			if (op.type == DIAG) {
				int diff = 0, pos = 0,count=0;
				for (auto itb = sseConf_->bondsArray.begin();
						itb != sseConf_->bondsArray.end(); itb++) {
					if (sseConf_->spin[(*itb)[0]] != sseConf_->spin[(*itb)[1]]) {
						int sign = (pos >= numOfSites_) ? -1 : 1;
						diff += sign;
						count++;
					}
					pos++;
				}
				if ((op.bond%2) && (op.bond<numOfSites_))
				*XPos1++=(double(diff)/(double(count)))/sseConf_->g;
				else
				*XPos1++=double(diff)/(double(count));
				*XPos2++=0.;
			}
		});
	// Corr
	for (int p = 0; p < (expRank - 1); p++) {
		for (int m = 0; m < expRank; m++) {
			XCorr_[p] += (3. / 2.)
					* (2. * Xp1_[m] * Xp1_[(m + p + 1) % expRank]
							+ Xp1_[m] * Xp2_[(m + p + 1) % expRank]
							+ Xp2_[m] * Xp1_[(m + p + 1) % expRank]);
		}
	}
// Draw Tau

	for (vector<double>::iterator it = tau_.begin(); it != tau_.end(); it++)
		*it = X_.draw();
	std::sort(tau_.begin(), tau_.end());
	vector<double>::iterator it_tau = tau_.begin();
	vector<double>::iterator it_XCorr = XCorr_.begin();
	for (; it_XCorr != XCorr_.end(); it_XCorr++, it_tau++)
		for (int k = 0; k < numOfOmega_; k++)
			gOmega_[k] += (*it_XCorr) * cos(2. * PI * (*it_tau) * double(k));
	LF_ << (1. / (pow(beta_, 2))) * gOmega_;
}

bool LodounFlourieMeasure::checkDone() const {
	return LF_.checkDone();
}
void LodounFlourieMeasure::dumpToFile(hdf5RW &hdf5File) const {
	LF_.dumpToFile(hdf5File, "LF");
}
void LodounFlourieMeasure::readFromFile(hdf5RW &hdf5File) {
	LF_.readFromFile(hdf5File, "LF");
}
int LodounFlourieMeasure::getNumOfMeasure() const {
	return LF_.getNumOfMeasure();
}

// Measurement Class

Measurements::Measurements() {

}

Measurements::~Measurements() {

}

void Measurements::init(SSEConf *sseConf, QmcParameters &qmcParameters,
		ofstream *outFile, MTGen *gen) {
// Init Measure
	outFile_ = outFile;
	allDone_ = false;
	sseConf_ = sseConf;
	ToMeasure toMeasure = qmcParameters.toMeasure.getVal();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	rate_ = qmcParameters.rate.getVal();
	progBar_.init(totalNumberOfMeasure, 0, qmcParameters.endTime.getVal(),
			qmcParameters.deltaDump.getVal(), outFile);
	if (toMeasure["WIND_ENERGY"])
		measures_.push_back(new WindingEnergyMeasure(sseConf_, qmcParameters));
	if (toMeasure["LF"])
		measures_.push_back(
				new LodounFlourieMeasure(sseConf_, qmcParameters, gen));
	if (toMeasure["HCORR"])
		measures_.push_back(new HcorrMeasure(sseConf_, qmcParameters, gen));

	//	if (toMeasure["STAGGERED"])
//		measures_.push_back(new StaggeredMeasure(sseConf_, qmcParameters));
}
void Measurements::DoMeasurement(bool &endTime, bool &dumpTime) {
	static int step = 0;
	endTime = dumpTime = false;
	step++;
	if (!(step % rate_)) {
		allDone_ = true;
		vector<int> numOfMeasure;
		for (MeasureVectorIter itm = measures_.begin(); itm != measures_.end();
				itm++) {
			bool done = itm->checkDone();
			allDone_ = allDone_ && done;
			if (!done)
				itm->DoMeasurement();
			int cNumOfMeasure = itm->getNumOfMeasure();
			numOfMeasure.push_back(cNumOfMeasure);
		}
		int minNumOfMeasure = *(std::min_element(numOfMeasure.begin(),
				numOfMeasure.end()));
		progBar_.update(minNumOfMeasure, endTime, dumpTime);
	}
}
void Measurements::updateMinNumOfMeasure(int totalNumberOfMeasure) {
	int minNumOfMeasure = totalNumberOfMeasure;
	for (MeasureVectorIter itm = measures_.begin(); itm != measures_.end();
			itm++) {
		int cNumOfMeasure = itm->getNumOfMeasure();
		minNumOfMeasure = std::min(minNumOfMeasure, cNumOfMeasure);
	}
	progBar_.updateMinNumOfMeasure(minNumOfMeasure);
}
void Measurements::dump(string outPath) const {
	string fileName = outPath + "/Measure.h5";
	hdf5RW hdf5File;
	hdf5File.initNew(fileName);
	for (auto it = measures_.begin(); it != measures_.end(); it++)
		it->dumpToFile(hdf5File);
}
void Measurements::read(string outPath) {
	string fileName = outPath + "/Measure.h5";
	hdf5RW hdf5File;
	hdf5File.initRead(fileName);
	for (auto it = measures_.begin(); it != measures_.end(); it++)
		it->readFromFile(hdf5File);
}
bool Measurements::checkNotDone() const {
	if (allDone_) {
		(*outFile_) << "Done Measurements" << endl;
		return false;
	} else {
		return true;
	}
}
void Measurements::updateRate(int const rate) {
	rate_ = rate;
}
int Measurements::getRate() const {
	return rate_;
}
// Progress Bar

void ProgressBar::init(int totalNumberOfMeasure, int minNumOfMeasure,
		double simTime, double deltaDump, ofstream *outFile) {
	startTime_ = boost::chrono::system_clock::now();
	simTime_ = simTime;
	totalNumberOfMeasure_ = totalNumberOfMeasure;
	outFile_ = outFile;
	deltaDump_ = deltaDump;
	this->updateMinNumOfMeasure(minNumOfMeasure);
}
void ProgressBar::update(int numOfMeasure, bool &endTime, bool &dumpFlag) {
	double p = (double(numOfMeasure) / double(totalNumberOfMeasure_));
	if (p > pc_) {
		pc_ += deltaDump_;
		dumpFlag = true;
		Hour diff = boost::chrono::system_clock::now() - startTime_;
		if (diff.count() > simTime_)
			endTime = true;
		(*outFile_) << "Finished " << p * 100 << " % , Time - " << diff.count()
				<< " Hours " << endl;
	}
}
void ProgressBar::updateMinNumOfMeasure(int minNumOfMeasure) {
	pc_ = ceil(
			(double(minNumOfMeasure) / double(totalNumberOfMeasure_))
					/ deltaDump_) * deltaDump_;
	(*outFile_) << "Finished Till Now" << pc_ * 100 << " % " << endl;
}

